Selection and validation of appropriate reference genes for RT–qPCR analysis of Nitraria sibirica under various abiotic stresses

Background Nitraria sibirica Pall. is a halophytic shrub with strong environmental adaptability that can survive in extremely saline-alkali and drought-impacted environments. Gene expression analysis aids in the exploration of the molecular mechanisms of plant responses to abiotic stresses. RT–qPCR is the most common technique for studying gene expression. Stable reference genes are a prerequisite for obtaining accurate target gene expression results in RT–qPCR analysis. Results In this study, a total of 10 candidate reference genes were selected from the transcriptome of N. sibirica, and their expression stability in leaves and roots under different treatment conditions (salt, alkali, drought, cold, heat and ABA) was evaluated with the geNorm, NormFinder, BestKeeper, comparative ΔCt and RefFinder programs. The results showed that the expression stability of the candidate reference genes was dependent on the tissue and experimental conditions tested. ACT7 combined with R3H, GAPDH, TUB or His were the most stable reference genes in the salt- or alkali-treated leaves, salt-treated roots and drought-treated roots, respectively; R3H and GAPDH were the most suitable combination for drought-treated leaves, heat-treated root samples and ABA-treated leaves; DIM1 and His maintained stable expression in roots under alkali stress; and TUB combined with R3H was stable in ABA-treated roots. TBCB and GAPDH exhibited stable expression in heat-treated leaves; TBCB, R3H, and ERF3A were stable in cold-treated leaves; and the three most stable reference genes for cold-treated roots were TBCB, ACT11 and DIM1. The reliability of the selected reference genes was further confirmed by evaluating the expression patterns of the NsP5CS gene under the six treatment conditions. Conclusion This study provides a theoretical reference for N. sibirica gene expression standardization and quantification under various abiotic stress conditions and will help to reveal the molecular mechanisms that confer stress tolerance to N. sibirica. Supplementary Information The online version contains supplementary material available at 10.1186/s12870-022-03988-w.

a series of complex mechanisms to resist or adapt to adverse environmental conditions. In addition, N. sibirica can also serve as an economic plant in saline-alkali areas due to its medicinal value and edible fruit [3]. In recent years, the physiological mechanisms of resistance to stress in N. sibirica have been studied, including the distribution of Na + and K + , osmotic adjustment and changes in antioxidase activity [4][5][6]. However, knowledge of the molecular mechanism of N. sibirica stress resistance is relatively insufficient.
Gene expression analysis is the foundation for revealing gene function and can help to reveal the molecular mechanisms involved in plant stress resistance. RT-qPCR (real-time quantitative polymerase chain reaction) is the most common technique used to analyse gene expression because it is highly sensitive, cost-effective, reproducible and specific. RT-qPCR has been widely used to investigate the expression patterns of genes in diverse organisms or under different environmental conditions [7][8][9]. However, the accuracy and reliability of RT-qPCR results are affected by many factors, such as the quality of the RNA, the efficiency of cDNA synthesis, and the efficiency of amplification [10,11]. To obtain accurate gene expression results, the data generated by RT-qPCR should be normalized with the use of appropriate reference genes [7,11]. Housekeeping genes that maintain the cytoskeleton or participate in basic biochemical metabolism, such as actin (ACT ), β-tubulin (TUB), elongation factor 1-α (EF-1α) and glyceraldehyde 3-phosphate dehydrogenase (GAPDH), are usually selected as internal reference genes [12][13][14]. The expression levels of reference genes should ideally be consistent in different tissues and under different environmental conditions [15,16]. However, numerous studies have found that the stability of reference gene expression is different in different species, tissues and experimental treatments [17,18]. For example, ACT was a stable reference gene in Platycladus orientalis [19] but unstable in teak (Tectona grandis L.f.) [20]. Ubiquitin conjugating enzyme (UBC) and GAPDH were the most stable reference genes in the flowers of Iris germanica L. '00246' and 'Elizabeth' , while the most stable reference genes in the flowers of Iris germanica L '2010200' were TUB and UBC [21]. In Glehnia littoralis, expressed protein 1 (EXP1) and serine/threonine-protein phosphatase PP2A (PP2A) were the most stable reference genes under salt stress, while cyclophilin 2 (CYP2) and α-TUB were the most stable genes under MeJA treatment [22]. In addition, it has been reported that some newly characterized genes may be more stable than traditional housekeeping genes. Under salt and drought stress, the expression of the new reference genes RG1, RG3 and RG5 in poplar was more stable than that of 18S ribosomal RNA (18S rRNA), Actin and ubiquitin (UBQ) [23]. The expression stability of unigene 14,912 was higher than that of 18S rRNA, GAPDH and His3 in annual ryegrass (Lolium multiflorum) under saline-alkali stress [24]. These studies show that there is not a universal reference gene for all experiments and all plant species [16,25]. It is necessary to systematically screen and identify internal reference genes before carrying out RT-qPCR experiments in different species, treatment conditions and different tissues. To date, there have been no reports on the screening and evaluation of appropriate reference genes in N. sibirica under abiotic stress.
The objective of the present study was to determine the optimal reference genes for RT-qPCR analysis in N. sibirica under various treatment conditions, such as salt, alkali, drought, heat and cold stresses and ABA treatments. To accomplish this objective, we screened N. sibirica transcriptome data for candidate reference genes, and five statistical algorithms (geNorm [26], NormFinder [10], BestKeeper [27], comparative ΔCt [28] and RefFinder [29]) were applied to evaluate their expression stability for normalization. In addition, the expression of Δ 1 -pyrroline-5-carboxylate synthetase (NsP5CS) was analysed to demonstrate the effectiveness of the selected reference genes. Our work will facilitate further studies of gene expression in N. sibirica and will accelerate the understanding of the molecular mechanisms underlying stress-induced responses in N. sibirica.

Primer specificity and amplification efficiency test of candidate reference genes
A total of 10 genes were selected from transcriptome data and used as candidate reference genes. qPCR primers were designed based on these sequences. All primers showed a signal and clear band of the expected size, and there was no primer-dimer formation observed on gel electrophoresis (Fig. S1), indicating the specificity of the primer pairs. The single peaks presented in the melting curve assays of each gene further verified the specificity of the primer pairs (Fig. S2). The amplification efficiency of the 10 candidate reference genes ranged from 89.67% (His) to 113.81% (ACT11), and the correlation coefficients (R 2 ) varied from 0.984 to 0.999 (Table 1). These results indicated that the primers of the 10 genes met the standards for qPCR and could be used in subsequent experiments. The details of the candidate reference genes, primer sequences, amplification lengths and efficiencies, and correlation coefficients are shown in Table 1.

Expression levels of candidate reference genes in leaves and roots
The quantification cycle (Cq) value reflects the mRNA transcript level. Reference genes with higher Cq values are considered to have lower expression abundance. Based on RT-qPCR, we determined the expression levels of 10 candidate reference genes in root and leaf tissues under salt, alkali, drought, cold, heat and ABA treatments. As shown in Fig. 1 (Fig. S4). Therefore, screening suitable reference genes for specific experimental treatments and tissues is needed.

Expression stability of candidate reference genes
Four specialized analytical tools, geNorm, NormFinder, BestKeeper and the comparative ΔCt method, were used to examine the expression stability of the 10 candidate reference genes. Subsequently, the RefFinder tool was employed to evaluate the expression stability of all these candidate reference genes and select the most suitable genes.

geNorm analysis
geNorm evaluated the stability of all 10 candidate reference genes using the M value (reference expression stability measure). The default threshold of M is 1.5. A lower M value indicates more stable gene expression [26]. As illustrated in Fig. 2, the M values of all 10 candidate reference genes in every treatment were lower than the default limit of 1.5, and all the candidate reference genes had different levels of stability under different treatments. ACT7 and R3H, which had the same M values, were the most stable genes for salt-treated leaves (SL), drought-treated leaves (DL) and drought-treated roots (DR); TUB showed good stability for salt-treated roots (SR), alkali-treated roots (AR), heat-treated roots (HR) and abscisic acid-treated roots (ABR); and GAPDH was one of the most stable genes in SR, alkali-treated leaves (AL), heat-treated leaves (HL) and abscisic acid-treated leaves (ABL). EF-1α and His were determined to be the least stable reference genes across treatments. geNorm software also provides information on the optimum number of reference genes to be used in an experiment based on the pairwise variation between ranked genes (Vn/n + 1), and a cut-off of 0.15 (V value) is usually applied to determine whether an additional reference gene needs to be added [26]. If Vn/n + 1 > 0.15, then n + 1 reference genes need to be used; otherwise, only n reference genes are needed. As shown in Fig

NormFinder analysis
NormFinder software directly evaluates the stability of reference genes according to the variance within and between groups, with lower values indicating higher stability [10]. The stability values of the candidate reference genes in each treatment are shown in Fig. 4 and Table S2. The stability rankings of the 10 reference genes under different experimental conditions were not completely consistent. For example, ACT7 was the most stable reference gene in the SL samples, while TUB was ranked first in the SR and ABR samples, and GAPDH was the top ranked gene for the AL, DL and HR samples. The stability ranking of the 10 candidate reference genes generated with NormFinder was slightly different from that produced by geNorm in most of the samples. For example, His and ACT11 were the top ranked reference genes for the DR and CR samples in the NormFinder analysis, while they ranked third and fourth, respectively, in the geNorm analysis.

BestKeeper analysis
BestKeeper evaluates the stability of gene expression by calculating the standard deviation (SD) and coefficient of variation (CV) of the Cq values. A smaller SD and CV indicate better stability of gene expression. If SD > 1, the gene was considered unstable [27]. The  and AR: alkali-treated leaves and roots, respectively; DL and DR: drought-treated leaves and roots, respectively; CL and CR: cold-treated leaves and roots, respectively; HL and HR: heat-treated leaves and roots, respectively; ABL and ABR: abscisic acid-treated leaves and roots, respectively analysis results are listed in Fig. 4 and Table S3. GAPDH was the most stable gene in the AL samples. His and TUB were ranked as the most stably expressed genes in the SR and DR samples, respectively, but exhibited the lowest stability in the HL and SL samples, respectively. DIM1 was the most stable reference gene in the SL and HL samples, and ACT11 was the most stable reference gene in the HR and ABR samples. For the AR, DL, CL and CR samples, the most stable gene was TBCB. EF-1α was ranked as the least stably expressed gene in most of the samples, including the AL, AR, DL, DR, CL, HL, HR and ABL samples. Overall, the stability ranking of the 10 candidate reference genes generated with BestKeeper was different from that of geNorm and NormFinder.

Comparative ΔCt analysis
This method evaluated gene expression stability by calculating the mean standard deviation (SD) value of each gene. Here, the smaller the value is, the higher the stability [28]. In SL, AL, AR, DR, HL, and ABR samples, ACT7 was one of the two best reference genes for gene normalization. For ABL, DL, CL, and HR samples, one of the two most stable reference genes was R3H.
In addition, the comparative ΔCt analysis also showed that EF-1α was the least stably expressed gene in multiple samples, such as AL, DR, CL and HR samples ( Fig. 4 and Table S4).

RefFinder analysis
To reduce the influence of the limitations of a single algorithm, comprehensive stability rankings of the candidate reference genes were determined with the Ref-Finder program (https:// www. heart cure. com. au/ reffi nder/). The ranking of genes was computed as the geometric mean, and a lower geometric, a higher stability []. As shown in Fig. 4 and Table S5, ACT7 combined with R3H, TUB, GAPDH or His were the two most stable reference genes in the SL, SR, AL, and DR samples; R3H and GAPDH were suggested to be the most suitable combination for DL, ABL, and HR samples; for AR samples, the two most stable reference genes were DIM1 and His; and TUB and R3H were the most stable reference genes for ABR samples. TBCB and GAPDH were the most stable reference genes for HL samples; TBCB, R3H, and ERF3A were found to be the three most stable reference genes in CL samples; and the three most stable reference genes for CR samples were TBCB, ACT11 and DIM1. EF-1α was the least stable reference gene in most samples.

Reference gene validation
To verify the results generated through the analyses described above, the expression pattern of NsP5CS was examined in SL, DR, AL, CR, HL and ABR samples. P5CS is a rate-limiting enzyme in proline biosynthesis and plays an important role in controlling plant stress tolerance [30,31]. The two most stable reference genes (alone and in combination) and the two most unstable reference genes based on the comprehensive ranking results for each sample set were used in the validation test. As shown in Fig. 5, the relative expression levels of NsP5CS normalized with the two most stable reference genes in combination showed significant changes in SL, DR, CR and ABR samples and different trends among different sample sets. In the SL and CR samples, NsP5CS was continuously induced, and its relative expression reached the highest level at 48 h, with values of 9.08 and 6.45, respectively. In DR samples, NsP5CS was rapidly induced and reached a maximum expression level at 3 h and then fluctuated at a lower level with time. NsP5CS was also rapidly induced and reached a maximum value at 3 h and remained at a higher value at 6 h for ABR samples. The expression patterns of NsP5CS normalized with the two most stable reference genes alone or in combination exhibited consistency in SL, DR, CR and ABR samples. However, the expression patterns of NsP5CS were significantly different when unstable reference genes were used. In the AL and HL samples, the expression level of NsP5CS varied less, but there were also similar expression patterns of NsP5CS normalized by the two most stable reference genes (alone or in combination) and were significantly different by the two unstable reference genes.

Discussion
RT-qPCR is one of the primary methods used to detect gene expression and thereby help to reveal the response mechanism of plants under different stresses [32]. The accuracy of gene expression data obtained by RT-qPCR analysis relies on the use of reference genes [33]. The use of a stable reference gene(s) for normalization ensures that the target gene expression data generated by RT-qPCR are reliable. In contrast, the use of an inappropriate internal reference gene will lead to inaccurate results [34].
In this study, a total of 10 reference genes were selected from the transcriptome data of N. sibirica to assess their expression stability in leaves and roots under salt, alkali, drought, cold, heat and ABA treatment. As shown in Table 1 and S1, the 10 primer pairs were specific, generating a single peak in the dissolution curve. The amplification efficiency of the 10 candidate genes ranged from 89.67 to 113.81%, and the R 2 was > 0.98, within the practical range [22,35,36]. Hence, the PCR conditions were acceptable. The mean Cq of the 10 candidate reference genes varied from 21.02 to 27.25 (Fig. 1). Considering that genes with a moderate expression level (a Cq value of 15 to 30, for instance) provide the most accurate normalization [37], the genes selected in this study met this requirement for use in normalization. Moreover, as a narrow distribution range tends to represent low variability, TBCB, R3H and DIM1 should be considered stable reference genes, with SD values of 0.87, 0.99 and 1.01, respectively. The Cq value of the genes is for the total sample pool. Considering that the reference genes have different changes in different tissues and experimental conditions (Fig. S3, Table S1), the stability of genes needs to be analysed separately for specific tissues and specific treatments, and more tools are needed. geNorm, NormFinder, BestKeeper and comparative ΔCt are algorithms that have been widely used for determining the stability of reference genes based on statistical analysis [16,18,24]. In the present study, the rankings generated by geNorm, NormFinder and comparative ΔCt were more similar to each other and different from Relative expression of NsP5CS at 0, 3, 6, 12, 24, and 48 h following stress treatment using the selected reference genes for normalization. The validated reference gene(s) used as normalization factors were the two most stable reference genes (alone or in combination) and the two most unstable reference genes in different treatment groups. Error bars show the standard error calculated from three biological replicates those obtained by BestKeeper. For example, the results of geNorm, NormFinder and comparative ΔCt showed that ACT7 and TUB were the most stable reference genes in SL and SR samples, while BestKeeper ranked third in stability. A similar difference between BestKeeper and other programs has also been reported in Nitraria tangutorum [13], Poa pratensis L. [38] and Maiyuanjinqiu [39], likely the result of the different algorithms that these programs employ [39,40]. RefFinder is a comprehensive evaluation program that uses data from geNorm, NormFinder, Best-Keeper and comparative ΔCt to calculate the geometric mean of candidate reference genes and generate a comprehensive stability ranking [41,42]. In this study, five programs (geNorm, NormFinder, BestKeeper, comparative ΔCt and RefFinder) were used to evaluate the stability of candidate reference gene expression in N. sibirica. Similar methods have been used in previous studies for many species, such as Eriobotrya japonica [42], Angelica decursiva [16] and wheat (Triticum aestivum L.) [43]. It has been reported that two or more stable reference genes are necessary to obtain accurate results [44,45]. The geNorm program could determine the optimal number of internal reference genes for normalization based on the cut-off of 0.15(Vn/n + 1) [26,46]. In this study, the V2/3 values for the leaves and roots under salt, alkali, drought, heat, and ABA treatments were below 0.15, indicating that two reference genes would be sufficient under these conditions. According to this rule, it is better to use three genes for cold-treated leaves and roots. However, the 0.15 threshold is not a strict restriction [47]. In addition, considering that the V2/3 values of CL and CR samples were slightly different from 0.15 (0.005 and 0.03), it is not necessary to use three reference genes instead of two as part of the reference gene validation.
Actin genes, including ACT2/7, ACT8, and ACT11, have been widely used as reference genes [38]. In this study, the expression stability of ACT7 was higher than that of ACT11 in most of the samples. Moreover, ACT7 was evaluated as one of the most stably expressed internal reference genes in multiple sample groups, such as SL, SR, AL, and DR samples. In a previous study, ACT7 was used as the reference gene for normalization of salt response gene expression in N. sibirica [2,48]. This study further demonstrated the reliability of ACT7 as an internal reference gene for N. sibirica. However, the expression stability of ACT11 was higher than that of ACT2/7 in soybean. This may be caused by species peculiarity. At the same time, this result also indicates that the stability of individual members of reference gene families can be diverse, and the stability of other members cannot be determined based on one member's stability. R3H was a novel reference gene identified from N. sibirica with unknown function and was stably expressed in the SL, DL, ABL, ABR, CL, and HR samples. Some new genes with greater expression stability than traditional housekeeping genes have also been identified in Lolium multiflorum [24] and poplar [23]. Therefore, it is not necessary to consider only traditional housekeeping genes when screening stable reference genes. In this study, EF-1α was the least stable in AL, DR, CL, HL and HR samples, while the expression of EF-1α was one of the most stable reference genes in Nitraria tangutorum [13]. This result indicates that the expression level of the same gene may not be constant across multiple species [15]. In addition, GAPDH was the best reference gene in the AL, DL, ABL, HL and HR samples of N. sibirica, while it exhibited poor stability in SL and AR. The same reference gene in the same species may respond differently to different stresses, as in the case in Kentucky bluegrass [38] and A. stolonifera [9]. In summary, it is usually necessary to select reference genes specifically according to species, tissues, and treatments.
Numerous studies have reported significant variations in the expression levels of target genes when unstable reference genes were used as internal controls [41,49]. To confirm the accuracy of our results, the expression pattern of NsP5CS was studied in SL, AL, DR, CR, HL, and ABR samples. The expression of NsP5CS in SL, DR, CR and ABR samples was significantly induced, and similar expression patterns of PuP5CS2, AtP5CS1 and SbP5CS1 under salinity, drought, cold and ABA treatment were reported in Populus ussuriensis [50], Arabidopsis thaliana [51] and sorghum [52]. The expression of NsP5CS in HL samples changed very little, which is consistent with the report that AtP5CS is not induced by high temperature [53]. MsP5CS was significantly induced by alkaline stress in alfalfa leaves [54], while the expression of NsP5CS in our AL samples changed little. This result may be due to the different species and reference genes used in the expression. In addition, as shown in Fig. 5, in each experimental treatment group, if stable reference genes were used to standardize the expression of NsP5CS, although there would be some small differences, the expression patterns were very similar. As previous studies have shown, the small differences can usually be corrected with more reference genes [55]. However, if unstable reference genes were used for normalization, the expression pattern of NsP5CS would be significantly deviated. The results of the validation test show that the reference genes screened in this study were reliable.
In summary, we conducted a systematic analysis to support the selection of stable reference genes for RT-qPCR analysis in leaves and roots of N. sibirica under six different treatments for the first time. ACT7 combined with R3H, GAPDH, TUB or His were the most stable reference genes in SL, SR, AL, DR samples, respectively; R3H and GAPDH were the most suitable combination for DL and HR; DIM1 and His maintained stable expression in AR sample; TUB combined with R3H was stable in ABR sample; TBCB and GAPDH exhibited stable expression in HL group; TBCB, R3H, and ERF3A were stable reference genes in CL group; and for CR sample, the three most stable reference genes were TBCB, ACT11 and DIM1. These results will improve the accuracy of target gene expression quantification, facilitate the identification of stress-responsive genes and help to reveal the molecular mechanisms conferring stress tolerance in N. sibirica.

Sample preparation and treatments
Seeds of N. sibirica were collected from Keluke beach saline-alkaline soil in the Qaidam Basin of Qinghai Province, China (37 • 10'-37 • 20' N, 96 • 49'-97 • 37' E), and they were identified by Zhang Huaxin, Institute of Ecological Protection and Restoration, Chinese Academy of Forestry (Voucher: Zhang Faqi, 20,170,823, Kunming Institute of Botany, Chinese Academy of Sciences, China). The seed collection followed relevant institutional, national, and international guidelines and legislation. The seeds were soaked in warm distilled water (40-50 °C) for 24 h to promote germination [56]. After surface sterilization using 10% (m/v) sodium hypochlorite for 15 min, the seeds were sown on Murashige and Skoog (1962; MS) medium containing 3.0% sucrose and 0.68% (m/v) agar. Then, the seeds were cultured at 26 °C, a 14-h photoperiod, and a photon flux density of 45 µmol m − 2 s − 1 . Three-week-old N. sibirica seedlings were transplanted into glass test tubes with half-strength Hoagland's solution that was renewed every week and grown under the same conditions. After 4 weeks, seedlings with good and consistent growth were randomly collected for different experimental treatments. For salt, alkali and drought stress, seedlings were treated with 400 mM NaCl, 120 mM NaHCO 3 and 20% PEG 6000 (w/v, polyethylene glycol), respectively. For heat and cold stress, the seedlings were subjected to temperatures of 45 °C and 4 °C, respectively, in a lighted incubator. For hormone treatments, the seedlings were treated with 100 µM ABA. The leaves and roots were harvested at 3, 6, 12, 24 and 48 h of treatment. At the same time, leaves and roots of untreated seedlings were also harvested as a blank control. The harvested samples were washed with purified water, immediately frozen in liquid nitrogen and stored at − 80 °C for further analysis. Notably, all the above treatments were performed in three biological replicates, and three seedlings were collected from each replicate.

Total RNA extraction and cDNA synthesis
RNA extraction was performed using the E.Z.N.A. ® Plant Kit (OMEGA Bio-Tek, Doraville, GA, USA) according to the product manual. The integrity of the extracted RNA was assessed by 1.5% agarose gel electrophoresis. Samples with a 28S/18S ratio of approximately 2.0 and no dispersion could be used for subsequent experiments. The quality and purity of the extracted RNA was determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). The absorbance ratios of the RNA samples at 260/280 and 260/230 nm ranged from 1.8 to 2.0 and 2.0 to 2.4, respectively, indicating suitable quality for subsequent research. According to the HiFiScript gDNA Removal RT MasterMix (CoWin Biosciences, Beijing, China) kit instructions, 1 µg of high-quality RNA extracted from each sample was reverse transcribed into first-strand cDNA. The obtained 20 µL cDNA products were then diluted to a volume of 60 µL for RT-qPCR.

Candidate reference gene selection and primer design
Candidate genes with stable expression were first selected according to the preliminary N. sibirica transcriptome data (NCBI with the ID GSE113246 and PRJNA804704); then, the candidate genes were further screened by combining gene annotation and public reports [13,45]. Finally, a total of 10 genes were selected as candidate reference genes. These candidate reference genes included 8 traditional reference genes (ACT7, ACT11, TUB, ERF3A, TBCB, GAPDH, EF-1α and His) and two novel genes, DIM1 and R3H. Primer Premier 5.0 software [57] was used to design the candidate reference gene primers. The parameters were set as follows: melting temperature (Tm), 55-65℃; GC content, 40-60%; primer length, 17-26 bp; and amplification product size, 150-300 bp [21,58]. All primers (Table 1) were synthesized by RuiBiotech company (Beijing, China). All primers were verified by electrophoresis on a 1.5% agarose gel.

RT-qPCR conditions and amplification efficiency test
RT-qPCR was performed in 96-well plates by using a LightCycler ® 480II Real-Time PCR System (Roche Molecular Systems, Germany). cDNA was amplified by using UItraSYBR Mixture (CoWin Biosciences, Beijing, China). Each reaction mixture contained 10 µl of 2×UItraSYBR Mixture, 0.8 µl of cDNA, 0.4 µl of each primer (10 µM) and 8.4 µl of ddH 2 O. The amplification conditions were as follows: an initial step of 95 °C for 10 min, followed by 40 cycles of 95 °C for 10 s, 55 °C for 30 s and 72 °C for 32 s. The final melting curve was produced by shifting the amplification temperature from 65 to 95 °C. RT-qPCR analysis of each sample was performed in triplicate, and template-free controls were included in parallel. The standard curve was constructed with a tenfold dilution (10, 100, 1000, 10,000) of a cDNA mixture comprising equal volumes of cDNA from all samples, and the amplification efficiency (E) and correlation coefficient (R 2 ) values of the primers were calculated using the standard curve. The E value of each primer pair was calculated by the curve slope using E = 10 (−1/slope) [59].

Stability analysis of candidate reference genes
To assess the stability of candidate reference genes, we first analysed the resulting RT-qPCR data using four software programs: geNorm, NormFinder, BestKeeper and comparative ΔCt method. Then, RefFinder (https:// www. heart cure. com. au/ reffi nder/) was used to generate a comprehensive ranking of the candidate reference genes according to data obtained by geNorm, NormFinder, BestKeeper and comparative ΔCt [42]. For geNorm and NormFinder analysis, the raw quantification cycle values (Cq) need to be converted into relative quantities, and for BestKeeper and comparative ΔCt algorithms, raw Cq values could meet the requirements.

Validation of reference gene stability
To validate the reliability of the selected reference genes, the relative expression level of NsP5CS was analysed in SL, AL, DR, CR, HL and ABR samples. The primer pair of NsP5CS (Table 1) was designed with Primer Premier 5.0 software. The RT-qPCR conditions were set up the same as the RT-qPCR conditions described above. The expression of NsP5CS was calculated using (E target ) ΔCq  [60] with the two worst and two best reference genes (alone and in combination) obtained by the comprehensive assessment used as references.